Experimental evolution of hybrid populations to identify Dobzhansky–Muller incompatibility loci

Abstract Epistatic interactions between loci that reduce fitness in interspecies hybrids, Dobzhansky–Muller incompatibilities (DMIs), contribute genetically to the inviability and infertility within hybrid populations. It remains a challenge, however, to identify the loci that contribute to DMIs as causes of reproductive isolation between species. Here, we assess through forward simulation the power of evolve‐and‐resequence (E&R) experimental evolution of hybrid populations to map DMI loci. We document conditions under which such a mapping strategy may be most feasible and demonstrate how mapping power is sensitive to biologically relevant parameters such as one‐way versus two‐way incompatibility type, selection strength, recombination rate, and dominance interactions. We also assess the influence of parameters under direct control of an experimenter, including duration of experimental evolution and number of replicate populations. We conclude that an E&R strategy for mapping DMI loci, and other cases of epistasis, can be a viable option under some circumstances for study systems with short generation times like Caenorhabditis nematodes.

formation and maintenance of distinct biological species (Coughlan & Matute, 2020).Dobzhansky-Muller incompatibilities (DMIs) are now appreciated as an important genetic basis for the evolution of such intrinsic post-zygotic isolation in speciation and yet the characterization of their location and identity in genomes remains challenging (Castillo & Barbash, 2017;Coyne & Orr, 1998;Presgraves, 2010b).
More than 80 years after Dobzhansky's pioneering work, the genes responsible for reproductive isolation remain incompletely characterized across organisms (Castillo & Barbash, 2017).However, recent technological developments and the increasing accessibility of genome-wide sequencing approaches have allowed for the isolation of specific DMI genes in laboratory settings (Blackman, 2016;Maheshwari & Barbash, 2011;Powell et al., 2020;Schumer et al., 2014).Such approaches include the introgression methods first pioneered by Dobzhansky, genome-wide association studies (GWAS), quantitative trait locus (QTL) mapping, clinal analysis, and linkage disequilibrium (LD) analyses (Barton & Hewitt, 1989;De La Torre et al., 2015;Fitzpatrick, 2013;Gompert & Buerkle, 2011;Payseur, 2010).LD mapping, in particular, has proved especially useful in determining not only the location but the number of loci participating in incompatible interactions (Schumer et al., 2020;Schumer & Brandvain, 2016), as demonstrated in swordtail fish by quantifying hundreds of DMI pairs in the naturally occurring hybrid zones of two species (Powell et al., 2020;Schumer et al., 2014Schumer et al., , 2015)).Specifically, they exploited patterns of correlation between allelic ancestry and linkage in the hybrid genomes in what is termed an ancestry disequilibrium scan: unlinked markers of the same parental origin that demonstrate linkage disequilibrium in hybrid populations indicate the presence of an incompatibility pair.
An analogous experimental approach to ancestry disequilibrium mapping for elucidating DMIs involves evolve-and-resequence (E&R) experiments (Kofler & Schlötterer, 2014;Long et al., 2015).E&R, as applied to mapping speciation genes (Pereira et al., 2021), combines multi-generation experimental evolution of hybrid populations with Pool-Seq, a modified next-generation high-throughput DNA sequencing approach that economically allows for the sequencing from a population pool of many individuals (Ferretti et al., 2013;Kofler et al., 2011;Kofler & Schlötterer, 2014;Turner et al., 2010Turner et al., , 2011)).While hypotheses on DMIs have rarely been tested directly in E&R studies, E&R has been applied to explore speciation and epistatic genetic interactions more generally (White et al., 2020).For instance, E&R using Escherichia coli to characterize the type of epistasis at play permitted the conclusion that negative epistasis contributes to intraspecific variation (Maharjan & Ferenci, 2013).The potential for E&R to provide a means to successfully map interspecies DMIs across genomes, however, remains to be fully explored.
Given the substantial resource and time investments required to carry out E&R, simulations provide a cost-effective means to guide long-term wet lab experiments that aim to characterize DMI loci (Schumer et al., 2020).One such simulation study assessed the power to map DMIs in yeast populations using linkage analysis, finding that large sample sizes and certain statistical methods performed best in successfully detecting DMI loci (Li et al., 2013).More recently, simulations explored the likelihood of hybrid speciation by elucidating how dominance, recombination, and selection influence DMI resolution (meaning that one or both alleles involved in an incompatibility pair get purged), though the power to localize DMI loci per se was not assessed (Blanckaert & Bank, 2018).Simulations of DMI inference from ancestry disequilibrium scans showed how ongoing hybridization of hybrid populations to parental populations can lead to biased co-ancestry in the genome (Schumer & Brandvain, 2016).
Here, we conduct forward-time computer simulations to explore the power of experimental evolution to detect and localize loci that contribute to epistatic genetic incompatibilities.By mimicking the Caenorhabditis nematode system that has partially isolated species pairs (Bundus et al., 2018;Cutter, 2018;Dey et al., 2014;Woodruff et al., 2010), we ground these simulations in a biological context with the potential to guide future E&R experiments while also assessing general properties of the potential for E&R mapping of DMIs to succeed.To date, Caenorhabditis has been used in experimental evolution across tens to hundreds of generations to study the evolution of sex, mutation rates, development, pathogen interactions, population genetic theory, thermal sensitivity, and other traits (Gray & Cutter, 2014;Teotónio et al., 2017), but not to map DMIs in hybrid populations.In contrast to co-ancestry disequilibrium scan approaches (Schumer & Brandvain, 2016), we aim to detect signatures at individual loci by taking advantage of replicated hybrid populations within an E&R framework.By exploring a range of parameters that define intrinsic properties of a system as well as features manipulable by an experimenter, we identify conditions most amenable to successful mapping of DMI loci with E&R experiments.

| Evolving hybrid populations in SLiM
To test the power of E&R experimental evolution studies to map DMIs, we simulated the evolution of hybrid populations containing DMIs using SLiM software (3.3.1)(Haller & Messer, 2019).We created artificial autosomal genomes containing: (1) a single twolocus DMI or (2) a suite of 10 randomly placed two-locus DMIs, then tracked their evolution in populations that mimic the hybridization of two species.Each replicate population contained 1000 diploid individuals, which stayed constant across generations.Over time, genotypes that created DMIs were selected against and purged from the populations, leading toward fixation of one or the other parental allele which we then attempted to map from patterns of parallel allele frequency changes across replicated hybrid populations.Each individual's genome started heterozygous at each locus, reflecting the logic that F 1 individuals that initialized the simulations would have inherited a distinguishable allele from each parental species (Figure 1).In the two-locus model, each locus of the DMI locus pair was linked to a separate autosome, with 10,000 evenly spaced neutrally evolving loci across the two autosomes as markers.No new mutations entered the populations over time.In the multi-locus model, the 20 loci that comprise 10 DMI pairs were randomly placed along two autosomes among an additional 10,000 evenly spaced neutrally evolving loci (approximately 1 marker per 2 kb).Because DMI locus positions in the multi-locus model were randomly assigned, they may overlap with neutrally evolving loci in some populations.While real-world systems might have >10,000 SNP markers, the LD between nearby markers over the modest timeframe of E&R would motivate pruning to a much smaller subset of diagnostic markers.
Each 20.8 Mb chromosome mimics the typical length of a Caenorhabditis chromosome (Teterina et al., 2020).Simulated genomes mimic the chromosome architecture of Caenorhabditis nematodes with a central region of low recombination and flanking arms of higher recombination (Rockman & Kruglyak, 2009), such that 50% of each chromosome is a low-recombination domain that contains 2500 neutral markers.Experimental evolution studies in Caenorhabditis, including C. elegans and C. remanei, commonly are run for tens to hundreds of generations (Gray & Cutter, 2014;Teotónio et al., 2017).Previous crosses between C. latens males and C. remanei females demonstrated X-autosome (causing hybrid male sterility), mito-autosomal, and autosomal-autosomal incompatibilities (the latter two responsible for inviability) (Bundus et al., 2018).The simulations, however, only consider autosomeautosome incompatibilities that are otherwise neutral, meaning that they do not individually confer a fitness advantage over alternative alleles.We also did not simulate explicitly the acquisition or analysis of short read sequence data, however, but assume that the simulated neutral marker loci represent distinguishable variants extracted from an E&R experiment.In this way, the simulation output may be used to guide the design of experiments with organisms such as C. latens and C. remanei to map hybrid incompatibilities between them.
We traced the evolutionary fate of alleles at DMI loci under a range of epistatic selection coefficients (s = 0.05, 0.2, 0.8 across all DMI pairs), interaction types (one-way or two-way interactions for the two-locus scenario), number of DMIs (multi-locus or twolocus scenarios with one-way interactions), and generation time sampling points (10-500 generations).Selection coefficients were chosen to allow us to investigate plausibly detectable effects ranging from mild to severe.For simulations of a single DMI pair, we encoded one DMI locus on an arm of one chromosome and the DMI partner locus on the center of the other chromosome (the highrecombination region "HRR" of 10 cM/Mb and low-recombination F I G U R E 1 (a) Schematic representation of the evolution of one-way and two-way DMIs from a single ancestral species.Our simulations were initiated with F1 hybrid individuals derived from the two parental species.(b) Simulated frequencies of alleles from parental species one over time in five replicate hybrid populations.Triangles and vertical dashed lines indicate positions of loci involved in a two-way DMI.While this example data was created using the same script as the simulated populations containing two-way DMIs in the results, it was generated to demonstrate the methods used and not for analysis.At Generation 1, all hybrid individuals in each population are F1 heterozygotes (allele frequency 0.5).In this example of two-way DMI, resolution of the DMI occurs upon fixation of alleles of both loci from parental species 1 or, alternately, upon fixation of alleles of both loci from parental species 2. The superimposed allele frequency profiles for the five replicate populations represented thus shown either extreme high or extreme low frequency near the DMI loci.
region "LRR" of 0.1 cM/Mb, respectively).For simulations of multiple DMI pairs, the DMI loci were randomly encoded throughout the two autosomes, irrespective of recombination regions; otherwise, the autosomes were constructed identically to those of the two-locus model.Parameter values were inspired by the biological range of what is known in Caenorhabditis species about chromosomal recombination landscapes and hybrid fitness (Bundus et al., 2018;Dey et al., 2014;Rockman & Kruglyak, 2009;Teterina et al., 2020).The simulated recombination rates are slightly more extreme than empirical averages, allowing these simulations to provide both high and low bounds with respect to less extreme recombination environments.
In modeling one-way incompatibilities, we followed the classical depiction of DMIs that occur when each parental species contributes only one allele to a derived-derived DMI pair (Figure 1).
Alternate allele combinations at the two loci would be compatible, because they would include ancestral alleles shared by the two species (Figure 1).In two-way incompatibilities, by contrast, each allele contributed by a parent at a DMI locus is incompatible with both of the alleles that could be contributed by the other parent at the partner DMI locus.Two-way DMIs thus result in two incompatibility types involving four derived alleles at two loci (Figure 1), which can occur when distinct mutations arose and fixed at both loci in each species separately after they diverged from their common ancestor, or from negatively epistatic ancestral-ancestral and ancestralderived allelic interactions (Figure 1).For the two-locus model with two-way DMIs, we simulated just two dominance scenarios.In one scenario, all four DMI alleles at each locus were dominant, resulting in the expression of both negative interactions starting in the first generation.Fitness effects from each incompatibility were multiplicative.For example, in the moderate selection case, for s = 0.106 for both DMIs present in an individual, fitness was reduced by 20% (80% of ancestral fitness) versus the 10.6% reduction (89.4% of ancestral fitness) when only one negative interaction was present.In the second scenario, all four DMI alleles were recessive, resulting in the expression of the DMIs only in later generations.Unlike in the first dominance scenario, both negative interactions could not be expressed at the same time in a given individual, resulting in a maximum fitness reduction of 10.6% when s = 0.106.

| Estimating DMI locations
Our simulations mimicked an E&R experiment by evaluating the power to detect DMIs from parallel allele frequency changes across replicate populations.We considered experimental evolution of 5, 10, and 20 replicate populations used to map DMIs for a given parameter set.Using an R script that processed allele frequencies to assess parallel evolution across a given set of replicate populations, we attempted to infer DMI locations on the simulated autosomes, iterating this procedure for each of 500 independent sets of population replicates to assess the precision and accuracy of mapping estimation (github.com/ Cutte rlab/ Gener ating -DMI-Estim ates; Figure 1).All simulation runs were carried out via gnu-parallel on the University of Toronto's Niagara/SciNet supercomputer (Ponce et al., 2019).
Estimates of the DMI locus locations were based on allele frequencies that were consistently more extreme than what would be expected from genetic drift alone.In some cases, this procedure inferred no DMI.To infer the location of DMI loci, we searched for overlapping genomic regions with extreme frequency values among the replicate populations.Allele frequencies were deemed extreme in a given replicate population when they exceeded values plausibly caused by genetic drift alone.The drift cutoff to define extreme allele frequencies was determined from 2.5th and 97.5th percentiles of 1000 neutral simulations that otherwise resembled the nonneutral simulations, including linkage structure, except that the alleles at the "DMI loci" exerted no fitness impacts.After counting the number of replicate populations with an extreme value at a given locus, we calculated with stats package in R the binomial probability of retrieving the observed number of extreme populations for each locus given the number of replicate populations considered (i.e., 5, 10, or 20 replicate populations), flagging loci with probabilities <.05 after false discovery rate correction.
In the two-locus scenario, genomic positions within the two-log ~95% confidence interval of the maximum negative log-binomial p-value for each chromosome were included to define the edges of the DMI location.If no loci significantly demonstrated an extreme number of populations with extreme allele frequencies for one or both chromosomes, then no DMI was estimated for that set of replicate populations (thus contributing to imperfect "detectability").We then assessed whether the estimated DMI locations were accurate based on whether they overlapped with the correct DMI positions input into the simulations, with accuracy defined as the proportion of correct detections regardless of the size of the detected region.The precision of the estimates was then calculated as the length of the regions that accurately overlapped true DMI loci divided by the length of the chromosome, subtracted from 1.This process was repeated 500 times for each distinct parameter set to assess the overall variability, detectability, accuracy, and precision of mapping DMI locations in the simulated genomes.We report the accuracy for a parameter set by tallying the number of estimates that overlapped the true DMI loci and dividing by the number of estimates made (usually 500, except when detectability was imperfect) and the precision for the parameter set by calculating the median relative length of the estimated DMI intervals.
In the multi-locus DMI simulations, we marked estimated DMI locus positions as blocks of adjacent loci with flagged probabilities <.05, provided that they included at least two loci.If two or more blocks of loci occurred within 6240 bp of each other (a span of 2.5 neutral loci), they were merged into a single estimate of DMI locus position.We report the number of loci inferred to

| Weak-to-moderate strength one-way DMIs can be mapped feasibly with E&R
When considering a single pair of DMI loci that express one-way incompatibility, the estimated map positions most accurately overlap the true DMI loci for intermediate durations of experimental evolution (Figure 2a).The precision of those map positions, however, increases with the number of generations of experimental evolution (Figure 2b).These general trends hold for weak to strong fitness effects, dominance regimes, and high-versus low-recombination environments; as expected, mapping precision is greater in highrecombination regions for a given duration of experimental evolution (Figure 2b).The dominance of interactions showed at most a negligible effect on overall mapping accuracy and precision (Figure 2).
Mapping accuracy and precision, however, are most sensitive to the duration of experimental evolution when selection is strong (s = 0.8), showing best accuracy for durations ≤100 generations and best precision for durations ≥100 generations.With weak selection (s = 0.05), by contrast, the relatively flat sensitivity to duration of experimental evolution yields the best trade-off between accuracy and precision between 40 and 200 generations.For both strong and weak selection on DMI loci, the precision of mapping is extremely poor and highly variable if populations experience 40 generations or fewer (the mapping interval often comprises >70% of the length of the chromosome), with the inference method unable to detect any DMI loci at all in 16% of runs after 10 generations with weak selection (Figure S2).
Intermediate selection on DMI loci (s = 0.2), however, permits robust accuracy and precision of DMI locus mapping from 10 to 100 generations in high-recombination regions, with mapping intervals spanning as little as 8% of the chromosome length (Figure 2b).Overall, an E&R approach may be most suitable to mapping pairwise one-way DMI loci of weak-to-moderate effect in a timeframe of 40-100 generations of experimental evolution.

| Accurate and precise mapping of most two-way DMIs with E&R requires 100 generations or more
When considering a single pair of DMI loci that express two-way incompatibilities, such that both alleles at one locus experience reciprocal negative interactions with both alleles at the other locus (Figure 1), we observed that both precision and accuracy of the estimated map positions tend to increase with the number of generations of experimental evolution (Figures 3 and 4).This pattern contrasts with the trends for one-way DMI loci that showed a trade-off between accuracy and precision (Figures 3a and 4a cf. Figure 2a).The improving accuracy and precision of two-way DMI mapping with duration of experimental evolution holds for weak to strong fitness effects, dominance regimes, and high-versus lowrecombination environments.Similar to one-way DMI loci, mapping precision is greater in high-recombination regions for dominant DMIs for a given duration of experimental evolution (Figure 3b).For instance, when selection is strong, we observed that mapping precision in high-recombination regions reaches a peak value of 96.5% versus a peak value of 83.1% for low-recombination regions.For recessive DMIs, however, a discernible beneficial influence of recombination on mapping precision holds true only when selection is strong (Figure 4b).More generally, we observe a clear interaction between dominance and strength of selection on the accuracy and precision of the estimated map positions for two-way DMI loci.
When selection is moderate or weak on incompatibilities, mapping precision tends to be better for dominant than for recessive DMI loci for any given duration of experimental evolution and, in abso-   2) to account for multiplicative fitness effects to retain a consistent maximum fitness reduction of 0.05, 0.2, or 0.8.There is no low-recombination point for the parameter set s = 0.025 at 80 generations because there were no accurate estimates retrieved for this parameter set (median accuracy = 0), and precision of mapped DMI regions are only calculated across accurate estimates.
but the shortest durations of experimental evolution and with only modestly reduced accuracy for durations of 40 generations and less.When recombination is high, we observed that mapping precision reaches a value of 96.1% by generation 100 for recessive DMIs (Figure 4b), compared to 78.4% for dominant DMIs (Figure 3b).

Strong selection on two-way recessive DMIs in high-recombination
regions thus provide the most promising circumstances for mapping two-way DMI loci, with nearly 90% precision and accuracy after just 20 generations of experimental evolution (Figure 4).However, most other circumstances that we explored require 100 generations or longer for high accuracy and precision in mapping two-way DMI loci.

| Complex trade-offs between true-and false-positive mapping with many pairs of DMI loci
Finally, we considered genomes with 20 DMI loci (10 pairs) that express one-way incompatibilities, such that just one allele at each locus experiences reciprocal negative interactions with one allele from another locus.We observed that the median size of estimated map locations tends to decrease with the number of elapsed generations (Figure 5a), consistent with the pattern of increasing precision observed in simulations of just a single pair of one-way DMI loci (Figure 2b).Median estimated size of DMI regions is not very sensitive to the range in number of replicate populations that we explored, though simulations of the fewest replicate populations As the regions spanning estimated DMI locations get narrower as generations elapse, the number of true DMI loci detected also increases (Figure 5b), reaching 90% or more of the 20 DMI loci after 100 generations with moderate selection (s = 0.02 for five replicate populations).However, the total number of estimated DMIs also increases with time, surpassing the true number of DMI loci, and leading to an accumulation of false-positives (Figure S2).Consequently, we observed that the fraction of true-positive estimates of DMI location exhibits a complex relationship with time and selection strength (Figure 5c).When selection is strong, true-positive detection of DMIs is highest in earlier generations when the regions span most of the genome and so have poor resolution.In later generations, true-positives decline to account for approximately 25% of mappings by generation 500 when the roughly 20 estimated DMI locations each represent about 0.1% of the length of a chromosome (Figure 5c).With strong selection, essentially all 20 DMI loci occur within the estimated regions, irrespective of the number of generations elapsed (Figure 5b).When selection is moderate (s = 0.02), the proportion of true-positive estimates for DMIs is highest with just 5 replicate populations until generation 150 (Figure 5c), as a result of the fewer estimates made with this lower powered design (Figure S2).However, the rapid growth of DMI estimates over time with moderate selection leads to a decline in true-positives down to ~14% by generation 500, irrespective of the number of population replicates (Figure 5c).As a result, early generations permit detection of only about ~45% of the 20 true DMI loci, reaching about ~100% of true DMI loci detected by generation 150 (Figure 5b).
With such weak selection, less than ~20% of true DMI loci (~4 of 20) are detected within the mapped regions at generation 10, reaching a maximum of 50% (10 of 20) by generation 500 and with only ~25% lying within estimated DMI regions even after this long timeframe (Figure 5b).

| DISCUSS ION
The genetic mapping of loci responsible for causing DMIs in interspecies hybrids represents an ongoing challenge (Maheshwari & Barbash, 2011;Presgraves, 2010b).The value of natural hybrid populations toward this goal, as in taxa like house mice and swordtail fish (Forejt et al., 2021;Moran et al., 2021;Powell et al., 2020), raises the possibility that experimental interspecies hybrid populations derived from experimental evolution might provide an underexploited technique to complement more traditional genetic mapping approaches (Gray & Cutter, 2014;Matute et al., 2020).We therefore simulated experimental evolution of hybrid populations to test for the power to detect DMI loci to inform the potential utility of an E&R experimental paradigm for this purpose.Overall, our observations indicate that mapping of DMIs with acceptable precision and accuracy using replicate hybrid populations is feasible, provided that the empirical system is capable of very rapid turnover of generations and that DMI fitness effects are neither too strong nor too weak.We found that the accuracy and precision of mapping trade-off in our simulations of one-way incompatibilities, such that an intermediate duration of experimental evolution is preferred for such circumstances.In contrast, mapping of two-way incompatibilities benefits from extension of the duration of experimental evolution for as long as possible.As expected, DMI loci in regions of high recombination always are more amenable to efficient mapping.We also observed that little benefit accrued for DMI detection with greater population replication over the range considered in our simulations (5-20 replicate populations) and that dominance interaction type exerted a strong influence just for two-way incompatibility detection.While an E&R approach to mapping DMI loci requires study systems with distinctive properties, with Caenorhabditis nematodes perhaps being especially suited to this technique (Castillo et al., 2015;Gray & Cutter, 2014;Teotónio et al., 2017), our findings support a place for experimental evolution as an additional tool in the toolkit of speciation genetics researchers (White et al., 2020).

| Weak-to-moderate effect DMIs are more feasible to map for one-way incompatibilities
We find that strength of selection, duration of experimental evolution, and recombination rate all affect the power to localize one-way DMI loci and show a trade-off between precision and accuracy that is especially sensitive to stronger fitness effects (Figure 2).Our results indicate that if selection in hybrids arises from just a small number of strong-effect DMIs, then mapping with an E&R experiment would require propagation for an infeasible number of generations in most common study systems, albeit possible for organisms like Caenorhabditis nematodes (i.e., the best accuracy-precision tradeoff of ~100 generations would take ~400 days).DMIs linked to high-recombination regions, of course, can be mapped with greater power for a given duration of experimental evolution (Blanckaert & Bank, 2018).A DMI locus pair that exerts weak to moderate effects, however, show more promise in yielding accurate and precise mapping resolution within 20-40 generations, with simulations showing the best trade-off at ~60 generations.Thus, when few DMI loci are anticipated, it may prove more expedient to pursue classic mapping strategies with controlled crosses.
The higher precision of map estimates for weak-and moderateeffect DMIs may be explained by the average time of DMI resolution: it takes longer to fix alleles that contribute to DMIs with weaker fitness effects (vertical dashed line in Figure S3).This slower resolution of DMI allows for more rounds of recombination to break up linked haplotype blocks surrounding the incompatibility loci.While this leads to estimates with higher precision, the resulting presence of fewer linked loci to draw information from compromises accuracy, demonstrating why the accuracies among weak selection coefficients do not exceed 90% in our simulations.
We did not find dominance relations between DMI loci to strongly affect the power to localize one-way DMI loci.In our design of an F 1 hybrid starting population, selection can quickly "see" both dominant and recessive incompatibilities starting in the F 2 generation.We expect dominance relations between DMI loci would exert a greater effect, however, when incompatibility alleles are not all present initially as heterozygotes.Incompatibility alleles starting at lower initial frequencies and in linkage disequilibrium might more closely match conditions simulated in other studies that report a stronger influence of DMI dominance interaction type (Blanckaert & Bank, 2018;Turner & Harr, 2014).

| Mapping of two-way incompatibilities benefits from long duration of experimental evolution
Two-way DMIs might be more common for species with greater overall divergence or for loci that experience rapid coevolution.In contrast to locus pairs involved in one-way incompatibilities, we observed both accuracy and precision of mapping two-way DMI loci to increase with the duration of experimental evolution (Figures 3 and   4).In two-way incompatibilities, all four alleles contribute to DMIs at the two loci and so selection affects both loci until both loci fix a single allele each.This requirement for fitness resolution in hybrid populations contrasts with one-way DMIs, which require fixation of only one locus, allowing polymorphic alleles at the partner DMI locus to then drift (Figure S3).The persistence of selection in the two-way scenario following fixation at one locus drives extreme allele frequencies at linked loci for both DMI loci to confer increased mapping power in later generations.Schumer et al. (2015) demonstrated a similar effect for what they referred to as symmetrical/ coevolved incompatibilities (two-way) versus asymmetrical/neutral incompatibilities (one-way) in their simulations of hybrid speciation.
Mapping of DMI loci that cause two-way incompatibilities also differed from one-way incompatibilities in that we observed a strong influence of dominance.We observed high precision and accuracy of mapping for dominant two-way DMIs with moderate-to-small fitness effects in high-recombination regions within a modest duration of experimental evolution (i.e., 20-40 generations, which could take ~40 to ~80 days in a real-world scenario applied to Caenorhabditis nematodes).Dominant two-way DMIs with the strongest fitness effects, however, suffer from poor precision until at least generation 150, even with high recombination.Recessive two-way DMIs, by contrast, show best mapping power when selection is strongto-moderate.For example, strong-effect recessive DMIs have both precision and accuracy values above 80% by generation 20 for high-recombination regions.These results indicate that mapping of two-way DMI loci may be most successful, for a given duration of experimental evolution, when the DMIs exert moderate fitness effects and lie in high-recombination regions.

| Abundant one-way DMI pairs can generate substantial false-positive mappings
While our analysis of single pairs of DMI loci proved instructive to assessing mapping power with an experimental evolution approach, the genomes of species are likely to encode many DMI loci (Masly & Presgraves, 2007;Schumer et al., 2014).Consequently, we assessed DMI mapping in the context of hybrid simulations involving 10 pairs of one-way incompatibilities.Reflecting the accuracy-precision trade-off that we documented for pairs of one-way incompatibilities, we observed that longer durations of experimental evolution did not enhance the incidence of true-positive mappings in the face of many When selection on DMIs is strong, this trade-off between accuracy and precision leads to the best compromise near generation 100, and essentially all true DMI loci are detectable.Small-tomoderate fitness effects of the DMIs, however, lead to many more false-positives in this timeframe and fewer true DMIs amenable to detection (i.e., lower accuracy), although the size of mapped regions are reasonably precise when they do overlap true DMI loci (mapping spans of 0.2% of chromosome length or less).These patterns indicate that if selection against many DMIs is suspected to be strong in real-world hybrid populations, then accurate and precise estimates of them could be achieved in a realistic amount of time in hybrid populations of organisms with rapid generation times.When selection against DMIs in hybrid genomes with many DMI pairs is weakto-moderate, mapping also can be feasible if one is willing to accept the presence of many false-positives.However, it is important to note that if selection against many DMIs is too strong, it could lead to population extinction.Population size was kept constant in these simulations and thus the potential for population extinction was not accounted for in these simulations, but could lead to lower power to detect DMIs in real-world settings.Similarly, our simulations excluded the potential for genetic drive elements or species-specific alleles that confer lab adaptive advantage, as well as fluctuations in population size, to influence allele frequency change.The presence of such factors could further obfuscate DMI mapping under real-world conditions.Moreover, hybrid genomes with high density of DMIs may lead to evolutionary reversion to a genomic state indistinguishable from just one of the parental species, as has been observed in experimental hybrid populations of Drosophila (Matute et al., 2020).This outcome could be exacerbated by higher genetic load in one of the parental species or by asymmetric potential for parental genotypes to be adaptive under lab conditions.Consequently, the power analyses of our simulations may overestimate mapping ability in the face of such additional real-world complications.
Finally, we found that the number of replicate populations exerted only a modest effect on the power to localize DMI loci in most cases with our approach, at least for the range of population replication considered in our simulations (5-20 populations).Consequently, excessive effort devoted to additional population replicates beyond 5 in real-world experiments often may not be worthwhile.This observation contrasts with the positive influence of replicates in traditional quantitative trait locus mapping in E&R experiments (Kofler & Schlötterer, 2014), and may reflect differences due to starting allele frequencies, additive versus epistatic selection, and alternative mapping inference techniques.

| Extensions to the investigation of DMIs with E&R
The populations modeled here simplify real-world DMIs, suggesting a number of potentially valuable extensions amenable to future work.For example, we only considered pairwise DMI interactions, but DMIs involving multiple interacting loci may be important in reproductive isolation (Dzur-Gejdosova et al., 2012;Satokangas et al., 2020;Schumer et al., 2014;Tiago et al., 2015;Turner & Harr, 2014).Moreover, it remains unresolved how often species ought to exhibit one-way versus two-way versus multi-way incompatibilities, though this may depend in part on the amount of accumulated divergence.Future work also could account for linkage between DMI loci, which is more likely when multiple genes are involved, rather than the encoding of partner loci on different chromosomes as in our simulations, or apply a co-ancestry disequilibrium scan approach to DMI detection in E&R studies (Schumer & Brandvain, 2016).Related to the issue of linkage is genome collinearity between species, such that the genomes of real species may show inversions and other rearrangements (Faria & Navarro, 2010;Fuller et al., 2018), as well as gene drives or alleles derived from one species that are adaptive under experimental evolution conditions, that could compromise the power to localize DMIs relative to simulations like ours that presume perfect synteny.Consequently, mapping of DMI loci with an E&R strategy should treat mapped regions as DMI candidates requiring further validation, as with any genetic mapping approach.Future work that explores the potential for E&R experiments to map DMIs will benefit from incorporating these issues, as well as other possible study designs (e.g., permitting gene flow from parental species), and direct contrasts of the efficacy of mapping DMIs via controlled crosses or QTL mapping.
Another important extension of this work would be sex-linkage or sex-limited expression of DMI loci (Coyne & Orr, 1998).Autosomal DMI loci with sex-limited expression will expose their fitness effects in only half of a population, so the power to detect them may resemble the small-effect DMIs simulated here.There is no real consensus, however, on what the fitness effects of DMI locus pairs should be, and thus it also will prove valuable to explore a broader spectrum of selection coefficients.Moreover, the DMI alleles that we modeled are neutral in their respective parental species' backgrounds, but adaptive alleles could also be modeled.Schumer et al. (2015) note that it is currently not known whether adaptive or neutral DMI alleles are more prevalent, and they modeled hybrid populations with both types of incompatibilities.In the adaptive case, alleles at distinct DMI loci are more likely to fix for the same parental species' allele combination (Schumer et al., 2015).Consequently, the outcomes of our two-way DMIs that also tend to fix DMI alleles in this way may capture some of the features of adaptive DMI allele introgression.
Another potential extension of this work would be to investigate whether the results could apply to detecting other forms of epistasis, as may be important among loci within species (Montooth et al., 2010;Whitlock et al., 1995).Detecting parallel allele frequency changes across replicate populations under specific parameter spaces might permit detection of loci implicated in meiotic drive and other intragenomic conflicts (Kofler et al., 2012;Rice, 2013).We welcome the future development of any alternative or more elegant approaches to detecting epistasis, whether within species or in the form of interspecies DMIs, in experimental populations that employ E&R study designs.

| CON CLUS ION
The ability to accurately and precisely detect DMI loci presents a challenge to characterizing the genetic underpinnings of reproductive isolation between species.Our exploration of the feasibility of experimental evolution with an E&R approach to map DMIs found that it will be most successful for modest fitness effects of DMIs caused by one-way incompatibilities, rather than two-way incompatibilities, that reside in high-recombination regions.DMI mapping in genomes that contain many locus pairs are subject to a trade-off between precision and accuracy that may be acceptable for strong-selection coefficients or if one is willing to accept the presence of many false-positives.Unlike the effects of differing selection coefficients and duration of experimental evolution, the number of replicate populations above five exerted only a modest effect on the power to localize DMI loci in most cases with our approach.Altogether, these findings demonstrate the plausibility of DMI mapping over feasible durations of experimental evolution with hybrid populations of organisms that have rapid turnover of generations, such as Caenorhabditis nematodes, to provide an addition to the toolkit of speciation genetics researchers.

For
simplicity in the multi-locus case involving 10 pairwise one-way DMIs, we simulated only one dominance scenario such that both DMI alleles of a pair acted dominantly.For one-way DMIs of the two-locus model, by contrast, we considered four dominance scenarios.First, both DMI alleles at each locus acted dominantly (dominant × dominant), resulting in expression of the incompatibility regardless of whether an individual was homozygous or heterozygous at a given DMI locus, and manifesting in the F 1 generation.Second, both DMI alleles were recessive (recessive × recessive), resulting in expression of the incompatibility only when an individual was homozygous for both DMI alleles and thus manifesting only in the F 2 and later generations.For the final two dominance scenarios (dominant × recessive), dominance differed by recombination region: the DMI allele in the high-recombination region was recessive and the DMI allele in the low-recombination region was dominant, or vice versa.The 10 DMI locus pairs conferred additive fitness effects, such that, for example, in the moderate selection case, each DMI decremented relative fitness by 0.02 for a maximum possible reduction of 20% across all 10 DMIs.
be mapped, the incidence of true-positive detection among the estimated DMI mappings, and the width of the regions spanning estimated DMI positions to assess the accuracy and precision of mapping DMI loci in genomes simulated to have 20 pairs of DMI loci.

F
I G U R E 2 (a) Accuracy and (b) precision of DMI mapping as a function of duration of experimental evolution in generations (log scale) for one-way incompatibilities.High-recombination region (HRR, red points) and low-recombination region (LRR, blue points) are shown for four different dominance landscapes (symbol shapes).Simulations map DMI loci based on parallel allele frequency changes among 10 replicate populations per parameter set.Accuracy is the number of mapped estimates that overlap true DMI loci divided by the total number of inferred estimates.Precision is one minus the median width (bp) of correct estimates divided by the length of one chromosome.Error bars around the median correspond to the 95% interquantile range (2.5th and 97.5th quantiles of mapped widths).Horizontal dashed lines correspond to accuracy and precision values of 50%, 80%, and 95%.
lute terms, mapping resolution of recessive incompatibilities is poor with fewer than 150 generations.Strong selection, by contrast, shows greater mapping precision for recessive two-way DMIs for all F I G U R E 3 (a) Accuracy and (b) precision values across generations of hybrid experimental evolution (log scale) for dominant × dominant two-way incompatibilities in the highrecombination region (red) and low-recombination region (blue).Accuracy, precision, error bars, and dashed lines as in Figure 2. Simulated selection parameters (s) per DMI are smaller for twoway DMIs relative to one-way DMIs (cf. Figure 2) to account for multiplicative fitness effects to retain a consistent maximum fitness reduction of 0.05, 0.2, or 0.8.F I G U R E 4 (a) Accuracy and (b) precision values across generations of hybrid experimental evolution (log scale) for recessive × recessive two-way incompatibilities in the highrecombination region (red) and low-recombination region (blue).Accuracy, precision, error bars, and dashed lines as in Figure 2. Simulated selection parameters (s) per DMI are smaller for twoway DMIs relative to one-way DMIs (cf.

(
five) tend to have widest genomic spans around estimated DMI locations.When selection on each DMI pair is weak (s = 0.005) or moderate (s = 0.02), mapping intervals are quite narrow (median estimated width does not exceed 0.2% of the genome).When selection is strong (s = 0.145), however, mapping resolution is entirely uninformative for the first 40 generations (Figure 5a).Specifically, we observed 100% of the genome typically to be implicated in incompatibility, regardless of the number of population replicates, though eventually reaching high precision after 60 generations or more (achieving mapping resolution for strong selection below 0.1% of chromosome length after 200 generations with 20 replicate populations; Figure 5a; Figure S1).

F
I G U R E 5 Power to detect DMI loci in experimentally evolving populations with 10 pairs of one-way DMI loci.(a) Median mapped estimate width (log scale fraction of genome), (b) proportion of true DMI loci detected, and (c) ratio of the number of true DMI loci mapped to the number of mapped regions, as a function of duration of experimental evolution (log scale generations) for 20 DMI loci that express one-way incompatibilities for 5, 10, and 20 replicate populations.Error bars indicate 2.5th to 97.5th quantiles.Square symbols in (c) indicate parameter sets for which mapped regions overlapped >1 true DMI on average.
DMI pairs in the simulated genomes.Allelic drift following fixation of one locus in a given DMI pair likely contributes to the inability of longer durations of experimental evolution to improve mapping confidence with many one-way DMI locus pairs encoded in genomes, based on our detection strategy that seeks to find parallel skews in allele frequency among replicate populations that exceed neutral expections.The median width of the region surrounding estimated DMI locations narrows with longer durations of experimental evolution and smaller selection coefficients.However, the number of estimated DMIs increases over time beyond the true number of DMI loci, leading the fraction of true-positives among estimated mapping locations to decline with the number of elapsed generations.